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We discuss the contributions to the Equation of State for the N£LO Skyrme pseudo-potential 
(£=2,3). We show that by adding 4th and 6th order gradient terms, it is possible to fairly reproduce 
the spin/isospin decomposition of an equation of state obtained from ab-initio methods. Moreover, 
by inspecting the partial-wave decomposition of the equation of state, we show for the first time a 
possible way to add explicit constraints on the sign of the tensor terms of the Skyrme interaction. 
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Introduction. A very successful nuclear physics model 
adopted to describe nuclear properties from drip-line 
to drip-line is the Nuclear Energy Density Functional 
(NEDF) theory [l|. Recently NEDF has been also ap¬ 
plied to study compact astrophysical objects as neutron 
stars M. To this purpose, the functional should be 
able to describe some infinite matter pseudo-observables. 
Among the different available functionals, the one de¬ 
rived from the effective non relativistic Skyrme interac¬ 
tion [ 5 ] is probably the most used by the nuclear struc¬ 
ture community. In its standard form @], the Skyrme 
functional consists of a linear combination of local den¬ 
sities up to second order in gradients plus a density 
dependent zero-range term which replaces the original 
three-body term [7|, |8j. Due to its very simple struc¬ 
ture, Brink and Vautherin Q performed the first nu¬ 
merical spherical Hartree-Fock (HF) calculations of even- 
even nuclei already in the early 70s. Once the struc¬ 
ture of the functional has been chosen, one should de¬ 
termine its coupling constants. This is done by building 
a penalty function and then performing a minimization 
in a multi-dimensional space [9j._Thc scientific collabo¬ 
ration named UNEDF-SciDAC 0.0 has investigated 
the role of the optimization procedure on the quality and 
predictive power of the standard Skyrme functional [ 13 - 
In their last article [13|, they have shown that there is 
no more room to improve the spectroscopic qualities of 
the standard Skyrme functional, by simply acting on the 
optimization procedure. Two possibilities are then avail¬ 
able: (i) following the DFT theory, where the primary 
building block is the functional that includes all corre¬ 
lation effects, (ii) following the self-consistent mean-field 
theory, where the major ingredient is an effective pseudo¬ 
potential and correlations are added afterwards using a 
beyond mean field approach. The main advantage of the 
second way is that we can follow a precise hierarchy in¬ 
cluding explicitly beyond-mean-field [if correlations to¬ 
wards the exact many-body ground state [lil ]. 


Following the second approach, Lesinski et al. have in¬ 
vestigated the inclusion of a tensor interaction [15t ] , while 
Sadoudi et al. considered an explicit central three-body 
term [§]. Carlsson et al. [l3] have studied the possibil¬ 
ity of including higher order derivative terms into the 
Skyrme functional. By means of the Density Matrix Ex¬ 
pansion method [17|, they have shown that it is possible 
to expand a finite range interaction in terms of a zero 
range like and such an expansion converges. In particular 
comparing the resulting binding energies of some selected 
nuclei, they have shown that using the standard Skyrme 
functional up to second order (N1LO), the average dis¬ 
crepancy is of the order of ss 10 MeV, while adding the 
4th order terms (N2LO) the error decreases by one or¬ 
der of magnitude. Despite these encouraging results, the 
major problem posed by this kind of approach lies in the 
strategy adopted to determine the higher order parame¬ 
ters. A badly constrained coupling constant can lead to 
numerical instabilities in finite nuclei calculations 0. 
In Ref. [19j , a numerical criterion based on the linear re¬ 
sponse (LR) formalism |2C} in symmetric nuclear matter 
(SNM) has been pres ented to avoid this kind of patholo¬ 
gies. In Refs. [21, 221, we have given explicit formulas for 
the N£LO cases discussed here, and we have also derived 
the necessary equations in coordinate space to be used 
in a numerical code. We also recall that the numerical 
code HOSPHE [2S| can already deal with higher order 
functionals (N3LO) to compute ground state energies of 
spherical even-even nuclei using the Harmonic Oscillator 
basis. Combining all these informations, it is thus possi¬ 
ble to build an optimization procedure to determine the 
parameters of the functionals avoiding by construction 
all bulk instabilities [24]. 


In this letter, we present a novel strategy to prop¬ 
erly constrain order by order the different parameters 
of the extended pseudo-potential by comparing with ab- 
initio calculations in infinite nuclear matter (INM). The 
use of INM pseudo-observables to constrain the Skyrme 
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functional has been already discussed by Chabanat et 
al. [2^|; in particular they included in the penalty func¬ 
tion the numerical values of the Equation of State (EoS) 
in pure neutron matter (PNM) of Wiringa et al. [26(. In 
recent years, we have seen a remarkable effort of the ab- 
initio community to improve the quality of their calcula¬ 
tions in INM 27-0]. All these data represent a valuable 
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opportunity to constrain the coupling constants of the 
N^LO functional. In fact, several ab-initio methods give 
us not only the general EoS in SNM, but also in each 
spin/isospin channel, and in some cases the partial wave 
contributions. In Ref. 
that a standard Skyrme functional is not flexible enough 
to reproduce these extra pseudo-data, otherwise it would 
be over-constrained and it would be incompatible with 
finite nuclei observables. There is in fact a strong dis¬ 
agreement, as shown in Fig. 6 of Ref. [0], not only in 
magnitude, but also in the sign of the energy per parti¬ 
cle in the different ( S,T ) channels as a function of the 
density. In the present letter, we show that the higher or¬ 
der derivative terms lead to a fair agreement with these 
pseudo-data. Moreover, by using the partial waves de¬ 
composition of the EoS [22|, it is possible to put con¬ 
straints on the sign of the tensor and spin-orbit N£LO 
terms. Such constraints are very useful since at present, 
because the existing parameterizations, based on ground- 
state observables of atomic nuclei, span in fact over a 
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wide range of magnitudes and signs [15|. 

Results. The extended Skyrme pseudo-potential 
(N£LO) contains central, spin-orbit and tensor com¬ 
ponents, whose general expressions have been de¬ 


duced in Refs. 21 


parameters t ^, Xq^ , 


l22j | . Besides the 
the parameters are 


zeroth-order 
labelled as 


4"4 4”4 4 "^ for the central part, and 4"\44 for 

the tensor part, where n = 2£ is the order of the 
N^LO expansion. With this notation, terms involving 
4 0 ) ,4 n) ,4 n) are even with respect to space exchange, 
and the remaining terms are odd. There is a single spin- 
orbit term, with coefficient Wn. which is originated from 
the second order expansion [30[. 

The EoS for the N£LO pseudo-potential in SNM can be 
easily obtained by means of a simple HF calculation [3Qj. 
Using the coupled spin-isospin basis (S', T) it is written 
as 
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where k F = (37 t 2 p/2) 1 / 3 , with p being the density, and 
y(S,T) potential energy per particle projected onto 

the different spin/isospin sub-spaces. The (S, T) channels 
in SNM up to 6-th order read 
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Neither the tensor nor the spin-orbit terms contribute to 
the different (S, T) channels. Using the expressions given 
in Eq. (HO, we have fitted the N2L0 and N3L0 param¬ 
eters on the Brueckner-Hartree-Fock (BHF) calculations 
of Baldo et al. [0] . The latter have been derived from 
the microscopic Argonne vlA nucleon-nucleon two-body 
interaction plus the Urbana model for the three-body 
term. For completeness, we have also considered the chi¬ 
ral effective field (%-EFT) calculations at low k of Hebeler 
et al. ini, which cover a narrower density interval than 
BHF results. 

We have found that no density-dependent term is 
needed to get a satisfactory fit. However, the result¬ 
ing parameters give a too low value for the effective 
mass, m*/m ~ 0.4. We have thus included, on top of 
the N£LO pseudo-potential, a density-dependent term. 
Actually, the effective mass does depend only on the 


combinations 34”^ + (5 


4x 


444 ”), but the presence 
of a density-dependent term produces a readjustment of 
the other parameters t\ , thus modifying the value of 
m*/m , an effect which has already been observed with 
the standard Skyrme interaction l32U34j . In general, such 
a density-dependent term should be obtained by apply¬ 
ing the method of 0,0,111 to three-, four-, ... body 
interactions. However, for the sake of simplicity, we 
have considered the standard Skyrme effective density- 
dependent term 4f 3 (l + XzP G )p a . Its contribution to all 
the formulae presented here can be immediately written, 
adding a similar contribution of 4°\ with the replace¬ 
ments 4 0) —> t 3 p a /6 and Xq—> X3. To keep a value 
of the effective mass m*/m ~ 0.7, we have fixed these 
parameters with two popular choices for a, namely 1/6 
and 1/3. It should be noticed that our conclusions do 
not depend on this particular choice. Furthermore, it is 
worth reminding that the present parameters should be 
considered only as a starting point for a more complete 
fit which should include finite nuclei observables. 

The N1L0 interaction can fit the ab-initio these re¬ 
sults only at low densities 0- On the contrary, a fair 
agreement is obtained with N2L0 and N3L0. In Fig. |T| 
are displayed as a function of the density the results for 
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FIG. 1: N2LO fits of V (s ’ T) (p) to the BHF (full dots) and 
X-EFT (full diamonds) results. 


the four (S, T) channels as obtained at the N2LO order. 
Dots and diamonds are the BHF 271 and y-EFT [36| 
results, respectively. Both are in reasonable agreement, 
except for the (1,1) channel. The N3LO fit of BHF chan¬ 
nels stay on top of the N2LO, pointing towards a conver¬ 
gence of the fits in that density range of values. This is 
not the case with the %-EFT results, because the data at 
our disposal are limited to density values <0.33 fm” 3 . 
The extrapolation to higher density values lead to very 
different N£LO, so that the fitted parameters are not reli¬ 
able enough. In the following we thus only refer to the fit 
to BHF results. In Tab. HI are given the numerical values 
of the central parameters of the pseudo-potential, here 
called VLyB6£ and VLyB3^, where l refers to the N£LO 
order, and the density dependent choice is a = 1/6,1/3, 
respectively. 


TABLE I: The parameters of the VLyB6£ and VLyB3£ 
pseudo-potentials, fitted using a standard density dependent 
term. Missing entries are zero, and n.d. stands for not deter¬ 
mined (see text). 



a = 1/6, tz 

=13763 MeV fm 3(1+a) , *3 = 1 

a = 1/3, £ 3 =8000 MeV fm 3(1+a ->, x 3 = 

1 


VLyB62 

VLyB63 

VLyB32 

VLyB33 

n i 

t\ n> [MeVfm 3+ "] 

(n) 

x\ 

tf ] [MeVfm 3-1 "" 

] x\ n> 

t\ n) [MeVfm 3+n ] 

x\ n> 

t\ n) [MeVfm 3+n ] 

(n) 

x\ 

0 0 

-2394.15 

0.632433 

-2500.49 

0.598658 

-1491.46 

0.409966 

-1660.94 

0.395793 

21 

-19.381 

35.182 

337.793 

-2.18967 

133.037 

-3.97968 

702.317 

-0.534133 

2 2 

513.2670 

-1.01914 

587.22 

-1.08109 

513.2670 

-1.01914 

587.22 

-1.08109 

41 

9.63577 

3.65615 

-116.489 

-0.477614 

-14.6617 

-0.745637 

-215.685 

0.201958 

4 2 

-65.3664 

-1.22006 

-108.342 

-1.33548 

-65.3664 

-1.22006 

-108.342 

-1.33548 

6 1 



2.84983 

-0.161801 



4.54218 

0.271071 

6 2 



1.23644 

-1.51102 



1.23644 

-1.51102 


A 71 ) 

i'O 

A™) 

t'e 

A n ) 

Lo 


A™) 
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A 71 ) 

t'e 

2 

401.816 

n.d. 

608.126 

n.d. 

401.816 

n.d 

608.126 

n.d. 

4 

-12.0604 

294.851 

-72.78 

-504.419 

-12.0604 

-225.429 

-72.78 

-221.001 

6 



20.7379 

107.776 



20.7379 

57.3209 


Wo=245.741 [MeVfm 3 ] 

W 0 =241.748 [MeVfm 3 ] 

W 0 =245.741 [MeVfm 3 ] 

W 0 =241.748 [MeVfm 3 ] 


Let us comment now on the fitted central parame¬ 
ters. A first trend is that within a given N£LO pseudo¬ 
potential, the parameters are always smaller, in 

absolute value, than t\ , which is a good signal about 
the convergence of the expansion in gradients. How¬ 
ever, the value of a given parameter can be very dif¬ 
ferent, even in sign, when going from N2LO to N3LO 
since the asymptotic behavior (i.e. the maximum power 
in the fit), is no longer governed by the same parame¬ 
ters. Thus, t j 4 -* in the N2LO fit, is constrained by the 
asymptotic behavior when we deal with N2LO, but it 
can take a very different value when higher order terms 
are considered in N3LO. By inspecting Eqs ([2|)-(pl|) one 


can see that the density dependent term only enters into 
the fit of channels (0,1) and (1,0), and consequently 
the space-odd parameters are independent of 

the choice of £*,£3,2:3. The N3LO space-even parame¬ 
ters are of the same sign and order of magnitude for 
the two selected density dependent terms. Using the 
central parameters given in Tab. HI we have calculated 
some basic SNM properties at the saturation density p sat . 
With VLyB62 interaction, we obtain the following values: 
p sat = 0.179 fm -3 , E/A = -15.50 MeV, m*/m = 0.789, 
and K a0 = 220.8 MeV. With VLyB63 these values are 
0.169 fm -3 , -15.72 MeV, 0.762, and 201.9 MeV, respec¬ 
tively. The interactions VLyB32 and VLyB33 give com- 
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parable results. 

To determine the spin-orbit and tensor contributions, 
we have to go one step further, and project the potential 
energy terms onto the (J,L,S,T) sub-spaces, where J 
and L are the total and orbital angular momenta, re¬ 
spectively. We have used the partial waves provided 
by the BHF microscopic calculations [37]. In Ref. [22] 
we have given explicit expressions for the potential en¬ 
ergy projections V( 2S+1 Lj), where we use the standard 
spectroscopic notation. The previous (5, T)-projections 
are obtained by summing up all these partial waves, and 
are independent of the spin-orbit and tensor parameters. 
This is also the case for the spin-singlet partial waves, 
which are thus fixed from the fitted central parameters. 
On the contrary, triplet spin partial waves do depend on 
the spin-orbit and tensor terms, leading to the removal 
of the degeneracy between the different partial waves for 
given J. 

The triplet partial waves relevant for our discussion are 
the following 


V( 3 Si) = ^4° ) (1 + :e o 0) )p+ -^ t i\ l + x i ) )P k F ( 6 ) 
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( 2 ) 

diagonal elements appear and the te ; term connects only 

states with L = 0 with states with L = 2. Consequently, 
( 2 ) 

r e ' cannot be determined with this type of calculation. 

By keeping the central parameters previously fitted to 
the (S, T)-channels, it is possible to determine the tensor 
and spin-orbit parameters ti 2,4,6 \ ti 4,6 \ Wo by fitting di¬ 
rectly to the BHF partial waves. However, the BHF cal¬ 
culation we use for the fit provides in fact two mixtures, 
namely 3 5i + 3 D±, and 3 P 2 + 3 F 2 . We have decided to 
fit the tensor parameters ti 4 '*, ti 6 ' 1 to the former of these 
combinations, and the spin-orbit Wo and the tensor 
parameters to the partial waves 3 Po and 3 Pi. The com¬ 
parison with the remaining BHF partial waves is in fact 
a prediction. In Fig. [2] are displayed the comparisons ob¬ 
tained at the N3LO level. We don’t display the N2LO 
results, because in contrast to the (S, T) channels, the fits 
are not so good, with notorious differences, in particular 
for the 1 S’o and 3 Si + 3 Di waves. 
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It is worth noticing that they do not depend on the tensor 
( 2 ) 

parameter te ■ Indeed, for the equation of state, only 


FIG. 2: (Color online) Comparing the partial waves 
V( 2S+1 Lj). The dots are the BHF results, and the lines are 
the results of our fit with VLyB3. 

One can see that the N3LO results are quite satisfac¬ 
tory, with the only appreciable difference for the 3 P2+ 3 E2 
coupled waves at high density values. Since the F wave 
receives contribution of order pkp only (see Eq. [TT]), it is 
clear that this discrepancy reflects a missing contribution 
of higher order (N4LO), as expected for a partial wave 
with £ = 3. 

Analogously to the central 4"\4"' 1 parameters, the 
spin-orbit Wo and tensor t^ parameters are fixed what¬ 
ever the choice made for the density-dependent term, as 
can be explicitly seen by inspecting Eqs m-m- The 
fits give practically the same Wo value (about twice the 
typical standard Skyrme values), with a positive sign. 
The same is true for whose N2LO and N3LO val¬ 
ues are of the same order of magnitude, being positive 
in both cases. Since this result is also in agreement with 
our previous work on Landau parameters [38j . we have 
now consistent evidence that the parameter to ; has to 
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be restricted to positive values. Finally, one can also see 
that within a given N£LO pseudo-potential, the 
parameters are always smaller, in absolute value, than 
£o?e, pointing again to a good convergence of the expan¬ 
sion in gradients. 

Conclusions. The inclusion of higher order derivative 
terms to the Skyrme pseudo-potential is of fundamen¬ 
tal importance to improve the spectroscopic qualities of 
the model, without loosing the previous expertise gained 
by the nuclear structure community since its very first 
applications in the early 70s. In this letter, we have 
demonstrated that the extra derivative terms can heal 
some major problems of standard Skyrme functionals. 
Performing a partial wave decomposition of the Skyrme 
pseudo-potential, we have been able to identify for the 
very first time, the contribution of the Skyrme tensor pa¬ 
rameters to the (S, T)-channels and partial waves. More¬ 
over, by fitting the BHF results, we have been able to 
derive a constraint on the sign of spin-orbit and tensor 
parameter to ' in total agreement with our previous work 
on Landau parameters. With the fitted parameters, we 
have calculated, as a function of the density, quantities 
as Landau parameters, and the EoS of polarized mat¬ 
ter, asymmetric matter, neutron matter and polarized 
neutron matter. No instabilities appear in the range of 
densities < 0.3 fm -3 which is relevant for finite nuclei. 
This set of parameters can thus be seen as a good start¬ 
ing point for a fit protocol which also includes constraints 
from finite nuclei. 

Since our analysis is mainly based on BHF calcula¬ 
tions, the question arises as to wether our method is of 
general reliability. We have seen that BHF and %-EFT 
results are in fair agreement at low densities, except in 
the (l,l)-channel. As y-EFT results are not yet available 
at high density values, we have not access to the asymp¬ 
totic behavior and a detailed comparison of the fits is 
not yet possible. However, our method enough to expect 
good fits to other ab-initio results, offering thus the pos¬ 
sibility of comparing among them, getting in touch with 
finite nuclei results through an effective NtLO pseudo¬ 
potential. 
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